# Code for Table I1: First Stage

rm(list = ls())

## ---------------------------------------
## Load Packages 
## ---------------------------------------
require('AER')
require('ivpack')
require('data.table')

## ---------------------------------------
## Load Data and Functions
## ---------------------------------------
load("../data0812.RData")

directory <- "../functions/"
functions <- list.files(directory)  
loadfunctions <- sapply(functions, FUN = function(x)source(paste0(directory, x)))

## ---------------------------------------
# Construct Instrument
## ---------------------------------------
data0812 <- constructIV(data0812)

## ---------------------------------------
#  Select Last Case Before Election
## ---------------------------------------
data0812 <- lastCase(data0812)

## ---------------------------------------
## Two Stage Least Squares Specifications
## ---------------------------------------

time.controls <- "as.factor(court_time1) + as.factor(court_time2) + as.factor(court_dow) + as.factor(court_shift) + as.factor(totOGS2)"
case.controls <-   "as.factor(any_drug_2) +  as.factor(any_violent_2) + as.factor(fire_arms_2) +  as.factor(any_rob_2) + as.factor(any_dui_2) + as.factor(prior_offender_2)"
case.controls_2 <- "as.factor(any_drug_2) +  as.factor(any_violent_2) + as.factor(fire_arms_2) +  as.factor(any_rob_2) + as.factor(any_dui_2)"
demo.controls <- "age_2012 + I(age_2012^2) + Female + as.factor(race) + vote2008 + as.factor(noteli08) + regis_before"
demo.controls_2 <- "age_2012 + I(age_2012^2) + Female + vote2008 + as.factor(noteli08) + regis_before"

outc.1 <- "vote2012"
endo.1 <- "pti"
inst.1 <- "judgeiv"

form.1 <- formula(paste(outc.1, "~", endo.1, "+" , time.controls, "|", inst.1, "+", time.controls))
form.2 <- formula(paste(outc.1, "~", endo.1, "+" , time.controls, "+", demo.controls, "|", inst.1, "+", time.controls, "+", demo.controls))
form.3 <- formula(paste(outc.1, "~", endo.1, "+" , time.controls, "+", demo.controls, "+", case.controls, "|", inst.1, "+", time.controls, "+", demo.controls, "+", case.controls))

form.1.fs <- formula(paste(endo.1, "~", inst.1, "+", time.controls))
form.2.fs <- formula(paste(endo.1, "~", inst.1, "+", time.controls, "+", demo.controls))
form.3.fs <- formula(paste(endo.1, "~", inst.1, "+", time.controls, "+", demo.controls, "+", case.controls))

## Prior case status
form.3_FO <- formula(paste(outc.1, "~", endo.1, "+" , time.controls, "+", demo.controls, "+", case.controls_2, "|", inst.1, "+", time.controls, "+", demo.controls, "+", case.controls_2))
form.3.fs_FO <- formula(paste(endo.1, "~", inst.1, "+", time.controls, "+", demo.controls, "+", case.controls_2))

## Race
form.2_B <- formula(paste(outc.1, "~", endo.1, "+" , time.controls, "+", demo.controls_2, "|", inst.1, "+", time.controls, "+", demo.controls_2))
form.3_B <- formula(paste(outc.1, "~", endo.1, "+" , time.controls, "+", demo.controls_2, "+", case.controls, "|", inst.1, "+", time.controls, "+", demo.controls_2, "+", case.controls))

form.2.fs_B <- formula(paste(endo.1, "~", inst.1, "+", time.controls, "+", demo.controls_2))
form.3.fs_B <- formula(paste(endo.1, "~", inst.1, "+", time.controls, "+", demo.controls_2, "+", case.controls))

## ---------------------------------------
## Row 1: Main Result - Full Sample
## ---------------------------------------

m1a1 <- ivreg(form.1, data = data0812)
m1a2 <- ivreg(form.2, data = data0812)
m1a3 <- ivreg(form.3, data = data0812)

m1d1 <- lm(form.1.fs, data = data0812)
m1d2 <- lm(form.2.fs, data = data0812)
m1d3 <- lm(form.3.fs, data = data0812)

## ---------------------------------------
## Row 2: Resource Deprivation
## ---------------------------------------

m2a1 <- ivreg(form.1, data = data0812[quants_inc == 1])
m2a2 <- ivreg(form.2, data = data0812[quants_inc == 1])
m2a3 <- ivreg(form.3, data = data0812[quants_inc == 1])

m2d1 <- lm(form.1.fs, data = data0812[quants_inc == 1])
m2d2 <- lm(form.2.fs, data = data0812[quants_inc == 1])
m2d3 <- lm(form.3.fs, data = data0812[quants_inc == 1])

m2b1 <- ivreg(form.1, data = data0812[quants_inc == 2])
m2b2 <- ivreg(form.2, data = data0812[quants_inc == 2])
m2b3 <- ivreg(form.3, data = data0812[quants_inc == 2])

m2e1 <- lm(form.1.fs, data = data0812[quants_inc == 2])
m2e2 <- lm(form.2.fs, data = data0812[quants_inc == 2])
m2e3 <- lm(form.3.fs, data = data0812[quants_inc == 2])

m2c1 <- ivreg(form.1, data = data0812[quants_inc == 3])
m2c2 <- ivreg(form.2, data = data0812[quants_inc == 3])
m2c3 <- ivreg(form.3, data = data0812[quants_inc == 3])

m2f1 <- lm(form.1.fs, data = data0812[quants_inc == 3])
m2f2 <- lm(form.2.fs, data = data0812[quants_inc == 3])
m2f3 <- lm(form.3.fs, data = data0812[quants_inc == 3])


## ---------------------------------------
## Row 3: Racially Disparate Impact
## ---------------------------------------
m5a1 <- ivreg(form.1, data = data0812[race == "Black"])
m5a2 <- ivreg(form.2_B, data = data0812[race == "Black"])
m5a3 <- ivreg(form.3_B, data = data0812[race == "Black"])

m5d1 <- lm(form.1.fs, data = data0812[race == "Black"])
m5d2 <- lm(form.2.fs_B, data = data0812[race == "Black"])
m5d3 <- lm(form.3.fs_B, data = data0812[race == "Black"])

m5b1 <- ivreg(form.1, data = data0812[race == "White"])
m5b2 <- ivreg(form.2_B, data = data0812[race == "White"])
m5b3 <- ivreg(form.3_B, data = data0812[race == "White"])

m5e1 <- lm(form.1.fs, data = data0812[race == "White"])
m5e2 <- lm(form.2.fs_B, data = data0812[race == "White"])
m5e3 <- lm(form.3.fs_B, data = data0812[race == "White"])

## ---------------------------------------
## Row 4: Prior case status
## ---------------------------------------
m4a1 <- ivreg(form.1, data = data0812[prior_offender_2 == 1])
m4a2 <- ivreg(form.2, data = data0812[prior_offender_2 == 1])
m4a3 <- ivreg(form.3_FO, data = data0812[prior_offender_2 == 1])

m4d1 <- lm(form.1.fs, data = data0812[prior_offender_2 == 1])
m4d2 <- lm(form.2.fs, data = data0812[prior_offender_2 == 1])
m4d3 <- lm(form.3.fs_FO, data = data0812[prior_offender_2 == 1])

m4b1 <- ivreg(form.1, data = data0812[prior_offender_2 == 0])
m4b2 <- ivreg(form.2, data = data0812[prior_offender_2 == 0])
m4b3 <- ivreg(form.3_FO, data = data0812[prior_offender_2 == 0])

m4e1 <- lm(form.1.fs, data = data0812[prior_offender_2 == 0])
m4e2 <- lm(form.2.fs, data = data0812[prior_offender_2 == 0])
m4e3 <- lm(form.3.fs_FO, data = data0812[prior_offender_2 == 0])

## ---------------------------------------
# First stage f-statistics

m1a1D <- summary(m1a1, vcov = sandwich, diagnostic = T)$diagnostics[1, 3]
m1a2D <- summary(m1a2, vcov = sandwich, diagnostic = T)$diagnostics[1, 3]
m1a3D <- summary(m1a3, vcov = sandwich, diagnostic = T)$diagnostics[1, 3]

m2a1D <- summary(m2a1, vcov = sandwich, diagnostic = T)$diagnostics[1, 3]
m2a2D <- summary(m2a2, vcov = sandwich, diagnostic = T)$diagnostics[1, 3]
m2a3D <- summary(m2a3, vcov = sandwich, diagnostic = T)$diagnostics[1, 3]

m2b1D <- summary(m2b1, vcov = sandwich, diagnostic = T)$diagnostics[1, 3]
m2b2D <- summary(m2b2, vcov = sandwich, diagnostic = T)$diagnostics[1, 3]
m2b3D <- summary(m2b3, vcov = sandwich, diagnostic = T)$diagnostics[1, 3]

m2c1D <- summary(m2c1, vcov = sandwich, diagnostic = T)$diagnostics[1, 3]
m2c2D <- summary(m2c2, vcov = sandwich, diagnostic = T)$diagnostics[1, 3]
m2c3D <- summary(m2c3, vcov = sandwich, diagnostic = T)$diagnostics[1, 3]

m4a1D <- summary(m4a1, vcov = sandwich, diagnostic = T)$diagnostics[1, 3]
m4a2D <- summary(m4a2, vcov = sandwich, diagnostic = T)$diagnostics[1, 3]
m4a3D <- summary(m4a3, vcov = sandwich, diagnostic = T)$diagnostics[1, 3]

m4b1D <- summary(m4b1, vcov = sandwich, diagnostic = T)$diagnostics[1, 3]
m4b2D <- summary(m4b2, vcov = sandwich, diagnostic = T)$diagnostics[1, 3]
m4b3D <- summary(m4b3, vcov = sandwich, diagnostic = T)$diagnostics[1, 3]

m5a1D <- summary(m5a1, vcov = sandwich, diagnostic = T)$diagnostics[1, 3]
m5a2D <- summary(m5a2, vcov = sandwich, diagnostic = T)$diagnostics[1, 3]
m5a3D <- summary(m5a3, vcov = sandwich, diagnostic = T)$diagnostics[1, 3]

m5b1D <- summary(m5b1, vcov = sandwich, diagnostic = T)$diagnostics[1, 3]
m5b2D <- summary(m5b2, vcov = sandwich, diagnostic = T)$diagnostics[1, 3]
m5b3D <- summary(m5b3, vcov = sandwich, diagnostic = T)$diagnostics[1, 3]

## ---------------------------------------
# Output

cont <- data.table(round(rbind(c(robust.se(m1d1)[2, c(1)], robust.se(m1d2)[2, c(1)], robust.se(m1d3)[2, c(1)]),
                               c(robust.se(m1d1)[2, c(2)], robust.se(m1d2)[2, c(2)], robust.se(m1d3)[2, c(2)]),
                               c(m1a1D, m1a2D, m1a3D),
                               c(m1a1$n, m1a2$n, m1a3$n),
                               rep(NA, 3),
                               c(robust.se(m2d1)[2, c(1)], robust.se(m2d2)[2, c(1)], robust.se(m2d3)[2, c(1)]),
                               c(robust.se(m2d1)[2, c(2)], robust.se(m2d2)[2, c(2)], robust.se(m2d3)[2, c(2)]),
                               c(m2a1D, m2a2D, m2a3D),
                               c(m2a1$n, m2a2$n, m2a3$n),
                               c(robust.se(m2e1)[2, c(1)], robust.se(m2e2)[2, c(1)], robust.se(m2e3)[2, c(1)]),
                               c(robust.se(m2e1)[2, c(2)], robust.se(m2e2)[2, c(2)], robust.se(m2e3)[2, c(2)]),
                               c(m2b1D, m2b2D, m2b3D),
                               c(m2b1$n, m2b2$n, m2b3$n),
                               c(robust.se(m2f1)[2, c(1)], robust.se(m2f2)[2, c(1)], robust.se(m2f3)[2, c(1)]),
                               c(robust.se(m2f1)[2, c(2)], robust.se(m2f2)[2, c(2)], robust.se(m2f3)[2, c(2)]),
                               c(m2c1D, m2c2D, m2c3D),
                               c(m2c1$n, m2c2$n, m2c3$n),
                               rep(NA, 3),
                               c(robust.se(m5d1)[2, c(1)], robust.se(m5d2)[2, c(1)], robust.se(m5d3)[2, c(1)]),
                               c(robust.se(m5d1)[2, c(2)], robust.se(m5d2)[2, c(2)], robust.se(m5d3)[2, c(2)]),
                               c(m5a1D, m5a2D, m5a3D),                               
                               c(m5a1$n, m5a2$n, m5a3$n),
                               c(robust.se(m5e1)[2, c(1)], robust.se(m5e2)[2, c(1)], robust.se(m5e3)[2, c(1)]),
                               c(robust.se(m5e1)[2, c(2)], robust.se(m5e2)[2, c(2)], robust.se(m5e3)[2, c(2)]),
                               c(m5b1D, m5b2D, m5b3D),                               
                               c(m5b1$n, m5b2$n, m5b3$n),
                               rep(NA, 3),
                               c(robust.se(m4d1)[2, c(1)], robust.se(m4d2)[2, c(1)], robust.se(m4d3)[2, c(1)]),
                               c(robust.se(m4d1)[2, c(2)], robust.se(m4d2)[2, c(2)], robust.se(m4d3)[2, c(2)]),
                               c(m4a1D, m4a2D, m4a3D),
                               c(m4a1$n, m4a2$n, m4a3$n),
                               c(robust.se(m4e1)[2, c(1)], robust.se(m4e2)[2, c(1)], robust.se(m4e3)[2, c(1)]),
                               c(robust.se(m4e1)[2, c(2)], robust.se(m4e2)[2, c(2)], robust.se(m4e3)[2, c(2)]),
                               c(m4b1D, m4b2D, m4b3D),                               
                               c(m4b1$n, m4b2$n, m4b2$n)
), 3))
# add column labels
colnames(cont) <- c("(1)", "(2)", "(3)")
# add row labels
cont$names <-  c("Main Effects",
                 "Std. Error",
                 "F-test", 
                 "Obs.",
                 "",
                 "Bottom Ter",
                 "Std. Error",
                 "F-test", 
                 "Obs.",
                 "Mid Ter",
                 "Std. Error",
                 "F-test", 
                 "Obs.",
                 "Top Ter",
                 "Std. Error",
                 "F-test", 
                 "Obs.",
                 "",
                 "Black",
                 "Std. Error",
                 "F-test", 
                 "Obs.",
                 "White",
                 "Std. Error",
                 "F-test", 
                 "Obs.",
                 "",
                 "Prior Case",
                 "Std. Error",
                 "F-test", 
                 "Obs.",
                 "No Prior Case",
                 "Std. Error",
                 "F-test", 
                 "Obs."
)

cont <- cont[, c("names", "(1)", "(2)", "(3)")]
cat("\nPrinting Table I1: First Stage...\n")
print(cont)
